#load the library
library(metafor)

#add the data for meta analysis
dat <- data.frame(
  study  = c("S3",  "S4a", "S4b", 
             "WP S1", "WP S2", "WP S3"),
  design = c("within", "between", "between", 
             "within", "between", "between"),
  
  m_low   = c(.61, 5.93, 6.05, .58, 3.45, 5.57),
  m_high  = c(.83, 5.30, 5.19, .94, 2.90, 5.16),
  
  sd_low  = c(0.67, 1.18, 1.10, 0.67, 1.59, 1.16),
  sd_high = c(0.79, 1.56, 1.51, 0.79, 1.41, 1.27),
  
  n_low   = c(100, 122, 80, 50, 112, 81),
  n_high  = c(100, 119, 80, 50, 101, 80))
  
##add cov for within-subject design here (NA for between-subjects)
cov = c(0.26, NA, NA, 0.18, NA, NA)
dat$r <- with(dat,
              ifelse(design == "within",
                     cov / (sd_low * sd_high),
                     NA))

##add the direction indicator
dat$higher_is_more_willing <- c(0, 1, 1, 0, 1, 1)
dat$direction <- ifelse(dat$higher_is_more_willing == 1, 1, -1)

#compute effect sizes
##between-subjects (S3, 4a, 4b, additional mediation, additional moderation) → SMD (Hedges g)
dat_between <- subset(dat, design == "between")

es_between <- escalc(
  measure = "SMD",
  m1i = m_high, sd1i = sd_high, n1i = n_high,
  m2i = m_low,  sd2i = sd_low,  n2i = n_low,
  data = dat_between
)

### Convert t-value from quantile regression to Hedges' g
compute_es5 <- function(t, df, n1, n2) {
  t_sign <- ifelse(t >= 0, 1, -1)
  r <- sqrt(t^2 / (t^2 + df)) * t_sign
  d <- 2 * r / sqrt(1 - r^2)
  es_5 <- escalc(
    measure = "SMD",
    di = d,
    n1i = n1,
    n2i = n2
  )
  return(es_5)
}

es_S5 <- compute_es5(t = -2.74, df = 186, n1 = 92, n2 = 98)
es_S5

##within-subjects (S3, pilot incentive compatible) → SMCR 
dat_within <- subset(dat, design == "within")

es_within <- escalc(
  measure = "SMCRH",
  m1i = m_high, sd1i = sd_high,
  m2i = m_low,  sd2i = sd_low,
  ni  = n_high,
  ri  = r,          # use the actual computed correlation
  data = dat_within
)

##combine all effect sizes
es_except_S5 <- rbind(es_between, es_within)
es_except_S5 <- es_all[match(dat$study, es_all$study), ]  # restore original order

es_S5_row <- data.frame(study = "S5", design = "between", m_low = NA, m_high = NA, sd_low = NA, sd_high = NA,
                        n_low = 92, n_high = 98, r = NA, higher_is_more_willing = 1, direction = 1,
                        yi = es_S5$yi, vi = es_S5$vi)
es_all <- rbind(
  es_except_S5[1:3, ],
  es_S5_row,
  es_except_S5[4:nrow(es_except_S5), ]
)
rownames(es_all) <- NULL

##adjust the direction of the effect size to be all consistent
es_all$yi_adj <- es_all$yi * es_all$direction
es_all$vi_adj <- es_all$vi
es_all

#run meta-analysis
res <- rma(yi_adj, vi_adj, data = es_all, method = "REML")
summary(res)

##forest plot
forest(res,
       slab = paste(es_all$study),
       xlab = "Effect size (high - low)",
       cex = 1.5)

